function Q0 = drawQ0(restr,phi,opt)
% Function draws Q Qdraws times from the space of orthonormal matrices 
% satisfying the equality restrictions represented in F and the sign 
% restrictions represented in signRestr. Note that this function assumes
% that the identified set is nonempty.
% Inputs are:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

F = restr.F;
f = restr.f;
S = restr.S;
s = restr.s;
Sigmatr = phi.Sigmatr;
Sigmatrinv = phi.Sigmatrinv;
signNorm = opt.signNorm;

N = size(Sigmatr,1); % Number of variables in VAR
mS = sum(s); % Number of sign restrictions

Q0 = zeros(N,N,opt.Qdraws); % Storage

for kk = 1:opt.Qdraws

    flag = 0;
    
    while flag == 0

        % Draw Q from space of orthonormal matrices satisfying equality 
        % restrictions.

        % Generate vector of independent standard normal random variables.
        z = randn(N,1); 

        qtilde = zeros(N);

        if isempty(F)

            qtilde(:,1) = z;

        else

            % Find rows of F restricting first column of Q (F1).
            F1 = F(1:f(1),:);

            % Compute residual from linear projection of z on F1'.
            [q,r] = qr(F1',0);
            qtilde(:,1) = z-F1'*(r \ (q'*z));

        end

        for jj = 2:N % For each column of Q

            if isempty(F)

                regMat = qtilde(:,1:jj-1);

            else

               % Find rows of F restricting jth column of Q (Fj).
               Fj = F((sum(f(1:jj-1))+1):sum(f(1:jj)),:);
               regMat = [Fj' qtilde(:,1:jj-1)];

            end

           % Generate vector of independent standard normal random variables.
           z = randn(N,1); 
           
           % Compute residual from linear projection of z on regMat.
           [q,r] = qr(regMat,0);
           qtilde(:,jj) = z - regMat*(r \ (q'*z));

        end

        Q = zeros(N);

        for jj = 1:N

            if signNorm == 0

                % Normalise diagonal elements of A0 to be positive.

                if Sigmatrinv(:,jj)'*qtilde(:,jj) == 0

                    % Randomly set sign of q_j to -1 or 1.
                    Q(:,jj) = (1-2*(rand > 0.5))*...
                        (qtilde(:,jj)./norm(qtilde(:,jj)));

                else

                    Q(:,jj) = sign(Sigmatrinv(:,jj)'*qtilde(:,jj))*...
                        (qtilde(:,jj)./norm(qtilde(:,jj)));

                end 

            else % Normalise diagonal elements of A0^(-1) to be positive

               if Sigmatr(jj,:)*qtilde(:,jj) == 0

                    % Randomly set sign of q_j to -1 or 1.
                    Q(:,jj) = (1-2*(rand > 0.5))*...
                        (qtilde(:,jj)./norm(qtilde(:,jj)));

               else

                   Q(:,jj) = sign(Sigmatr(jj,:)*qtilde(:,jj))*...
                       (qtilde(:,jj)./norm(qtilde(:,jj)));

               end

            end

        end

        % Check whether proposed draw satisfies sign restrictions.

        % Check whether proposed draw satisfies sign restrictions (if any).

        if mS > 0  % If sign restrictions

            signCheck = zeros(mS,1);
            restrCount = 0;

            for ii = 1:N

                for jj = 1:s(ii)

                    restrCount = restrCount + 1;

                    signCheck(restrCount) = S(restrCount,:)*Q(:,ii);

                end

            end

            % If sign restrictions satisfied, terminate while loop and 
            % save value of Q.
            flag = (min(signCheck) > 0);
            Q0(:,:,kk) = Q;

        else % If no sign restrictions

            Q0(:,:,kk) = Q;
            flag = 1;

        end

    end
    
end

end